% Matlab Code for replication of "Optimal Debt Maturity and Firm Investment"
% by J. Jungherr and I. Schott, September 2020
% Stationary distribution of firms

% For each Z level, create a fine DEBT grid
lenb2       = 500; 
debtgrids   = zeros(lenb2,lenz);
for zz = 1:lenz
    debtgrids(:,zz)   = linspace(PARMS.DEBT{zz}(1)*0.95,PARMS.DEBT{zz}(end),lenb2);
end

% Create interpolatable policy functions
ltdF        = cell(1,1);
defF        = cell(1,1);
for zz = 1:lenz
    ltdF{zz}    = griddedInterpolant(PARMS.DEBT{zz},ltd_pol(:,zz),interpmethod,extrapmethod);
    defF{zz}    = griddedInterpolant(PARMS.DEBT{zz},def_pol(:,zz),interpmethod,extrapmethod);
end

bp_pol_stay     = zeros(lenb2,lenz);
bp_pol_up       = zeros(lenb2,lenz);
def_pol2        = zeros(lenb2,lenz);
ltd_pol2        = zeros(lenb2,lenz);
err1            = zeros(lenb2,lenz);
err2            = zeros(lenb2,lenz);

for zz = 1:lenz
    for bb = 1:lenb2
        
        % Determine default probability at this point
        def_pol2(bb,zz)     = defF{zz}(debtgrids(bb,zz));
        
        bchoice             = ltdF{zz}(debtgrids(bb,zz));
        ltd_pol2(bb,zz)     = bchoice;
        
        % We want to know which grid point the firm will be on tomorrow
        [val, idx]          = min(abs(bchoice*(1-PARMS.GAMMA) - debtgrids(:,zz)));
        bp_pol_stay(bb,zz)  = idx;
        err1(bb,zz)         = val/(bchoice*(1-PARMS.GAMMA));
        
        [val, idx]          = min(abs(bchoice*(1-PARMS.GAMMA) - debtgrids(:,min(lenz,zz+1))));
        bp_pol_up(bb,zz)    = idx;
        err2(bb,zz)         = val/(bchoice*(1-PARMS.GAMMA));
        
    end
end


%% Transition matrix

n           = lenz*(lenb2*lenz); % number of possible combinations of states today and tomorrow
ntriplets   = n/lenz; % size of trans (in each dimension)
I           = zeros(n,1); % Row index
J           = zeros(n,1); % Column index
X           = zeros(n,1); % Value

ztrans                          = zeros(lenz);
ztrans(1:lenz+1:lenz^2)         = PARMS.rhoz;  
ztrans(1+lenz:lenz+1:lenz^2)    = 1-PARMS.rhoz;  
ztrans(end,end)                 = 1;


for bb = 1:lenb2  
    for zz = 1:lenz  
        
        today       = (zz-1)*lenb2 + bb;
        idx         = ((1:lenz)-1)*lenb2*lenz + today;
        I(idx)      = today; % row index
        bchoice     = ones(1,lenz)*bp_pol_stay(bb,zz);
        bchoice(min(lenz,zz+1)) = bp_pol_up(bb,zz);
        tom         = ((1:lenz)-1)*lenb2 + bchoice; % states tomorrow
        J(idx) = tom; % column index
        X(idx) = ztrans(zz,:)*(1-def_pol2(bb,zz))*(1-PARMS.x); % value
        
    end
end

trans   = sparse(I,J,X,ntriplets,ntriplets);
clear I J X
clear ztrans

%% Stationary distribution

M                   = 1; % normalization
entrants            = sparse(lenb2*lenz,1);
idx                 = (1-1)*lenb2 + 1;
entrants(idx)       = 1;
I                   = speye(lenb2*lenz); % sparse identity matrix
probst              = (I-trans')^(-1)*(entrants)*M;   
dist                = full(reshape(probst,lenb2,lenz));


%% Put policies on finer grid

cap_pol2    = zeros(lenb2,lenz);
lab_pol2    = zeros(lenb2,lenz);
epsbar_pol2 = zeros(lenb2,lenz);
std_pol2    = zeros(lenb2,lenz);
ktilde_pol2 = zeros(lenb2,lenz);
lev_pol2    = zeros(lenb2,lenz);
p_pol2      = zeros(lenb2,lenz);
ps_pol2     = zeros(lenb2,lenz);
mat_pol2    = zeros(lenb2,lenz);
eq_pol2     = zeros(lenb2,lenz);

for zz = 1:PARMS.lenz
    capF            = griddedInterpolant(PARMS.DEBT{zz},cap_pol(:,zz),interpmethod,extrapmethod);
    labF            = griddedInterpolant(PARMS.DEBT{zz},lab_pol(:,zz),interpmethod,extrapmethod);
    stdF            = griddedInterpolant(PARMS.DEBT{zz},std_pol(:,zz),interpmethod,extrapmethod);
    ktldF           = griddedInterpolant(PARMS.DEBT{zz},ktilde_pol(:,zz),interpmethod,extrapmethod);
    levF            = griddedInterpolant(PARMS.DEBT{zz},lev_pol(:,zz),interpmethod,extrapmethod);
    plF             = griddedInterpolant(PARMS.DEBT{zz},p_pol(:,zz),interpmethod,extrapmethod);
    psF             = griddedInterpolant(PARMS.DEBT{zz},ps_pol(:,zz),interpmethod,extrapmethod);
    matF            = griddedInterpolant(PARMS.DEBT{zz},mat_pol(:,zz),interpmethod,extrapmethod);
    epsbarF         = griddedInterpolant(PARMS.DEBT{zz},epsbar_pol(:,zz),interpmethod,extrapmethod);
    eqF             = griddedInterpolant(PARMS.DEBT{zz},eq_pol(:,zz),interpmethod,extrapmethod);

    % Apply to large grid
    cap_pol2(:,zz)  = max(capF(debtgrids(:,zz)),0);
    lab_pol2(:,zz)  = max(labF(debtgrids(:,zz)),0);
    epsbar_pol2(:,zz)  = epsbarF(debtgrids(:,zz));
    std_pol2(:,zz)  = stdF(debtgrids(:,zz));
    ktilde_pol2(:,zz)   = ktldF(debtgrids(:,zz));
    lev_pol2(:,zz)  = levF(debtgrids(:,zz));
    p_pol2(:,zz)    = plF(debtgrids(:,zz));
    ps_pol2(:,zz)   = psF(debtgrids(:,zz));
    mat_pol2(:,zz)  = matF(debtgrids(:,zz));
    eq_pol2(:,zz)   = eqF(debtgrids(:,zz));

end

